home *** CD-ROM | disk | FTP | other *** search
/ Best of www.BestZips.com (Collector's Edition) / Best of WWW.BESTZIPS.COM Collector's Edition (JCSM Shareware) (JCS Marketing).ISO / prgtools / tn2.zip / MOONUP.T < prev    next >
Text File  |  1996-11-15  |  11KB  |  498 lines
  1. %
  2. % "moonup.t" computes the times of moonrise and moonset 
  3. % anywhere in the world.
  4. %
  5. % Adapted from a BASIC program from SKY & TELESCOPE, July, 1989, page 78
  6. %
  7. % Enter north latitudes positive, west longitudes negative.  
  8. %
  9. % For the time zone, enter the number of hours west of Greenwich 
  10. % (e.g., 5 for EST, 4 for EDT).  
  11. %
  12. %   Sample program for the T Interpreter by:
  13. %
  14. %   Stephen R. Schmitt
  15. %   962 Depot Road
  16. %   Boxborough, MA 01719
  17. %
  18.  
  19. const PI : real := 2 * arcsin( 1 )
  20. const DR : real := PI / 180
  21. const K1 : real := 15 * DR * 1.0027379
  22.  
  23. var Moonrise, Moonset : boolean := false
  24.  
  25. type rvector : array[3] of real
  26.  
  27. program
  28.  
  29.     var k, zone : int
  30.     var jd, lat, lon, ph, t0, tz : real
  31.     var ra0, dec0, plx0 : real
  32.     var mp : array[3, 3] of real
  33.     var ra, dec, v : rvector
  34.     
  35.     prompt "Latitude, +N, -S (deg):"
  36.     get lat
  37.     prompt "Longitude, +E, -W (deg):"
  38.     get lon
  39.     prompt "Time zone +West of GMT, -East (hrs):"
  40.     get zone
  41.  
  42.     jd := julian_day - 2451545  % Julian day relative to Jan 1.5, 2000
  43.     show_input( lat, lon, zone )
  44.     
  45.     lon := lon / 360
  46.     tz := zone / 24
  47.  
  48.     t0 := lst( lon, jd, tz )
  49.     jd := jd + tz
  50.  
  51.     for k := 0...2 do
  52.  
  53.         moon( ra0, dec0, plx0, jd )
  54.         mp[k,0] := ra0
  55.         mp[k,1] := dec0
  56.         mp[k,2] := plx0
  57.         jd := jd + 0.5        
  58.         
  59.     end for
  60.  
  61.     if mp[1,0] <= mp[0,0] then
  62.         mp[1,0] := mp[1,0] + 2 * PI
  63.     end if
  64.  
  65.     if mp[2,0] <= mp[1,0] then
  66.         mp[2,0] := mp[2,0] + 2 * PI
  67.     end if
  68.  
  69.     ra[0]  := mp[0,0]
  70.     dec[0] := mp[0,1]
  71.  
  72.     put ""
  73.  
  74.     for k := 0...23 do
  75.  
  76.         ph := ( k + 1 ) / 24
  77.         
  78.         ra[2]  := interpolate( mp[0,0], mp[1,0], mp[2,0], ph )
  79.         dec[2] := interpolate( mp[0,1], mp[1,1], mp[2,1], ph )
  80.         
  81.         test( k, zone, t0, lat, mp[1,2], ra, dec, v )
  82.  
  83.         ra[0]  := ra[2]
  84.         dec[0] := dec[2]
  85.     
  86.     end for
  87.  
  88.     if ( not Moonrise ) and ( not Moonset ) then
  89.     
  90.         if v[2] < 0 then
  91.             put "Moon down all day"
  92.         else
  93.             put "Moon up all day"
  94.         end if
  95.  
  96.     else
  97.  
  98.         if not Moonrise then
  99.             put "No moonrise this date"
  100.         elsif not Moonset then
  101.             put "No moonset this date"
  102.         end if
  103.  
  104.     end if
  105.     
  106. end program
  107.  
  108. %
  109. % display latitude, longitude and time zone
  110. %
  111. procedure show_input( lat, lon : real, zone : int )
  112.  
  113.     var deg, min : int
  114.     
  115.     if sgn( zone ) = sgn( lon ) and zone ~= 0 then
  116.         put "WARNING: time zone and longitude are incompatible"
  117.     end if
  118.  
  119.     if abs( zone ) > 12 then
  120.         put "WARNING: invalid time zone"
  121.     end if        
  122.  
  123.     if abs( lat ) > 90 then
  124.         put "WARNING: invalid latitude"
  125.     end if
  126.  
  127.     if abs( lon ) > 180 then
  128.         put "WARNING: invalid longitude"
  129.     end if
  130.  
  131.     if lat > 0.0 then
  132.         deg := floor( lat )
  133.         min := floor( 60 * ( lat - deg ) )
  134.         put deg:4, ":",min:2, " N "...
  135.     else
  136.         deg := floor( -lat )
  137.         min := floor( 60 * ( -lat - deg ) )
  138.         put deg:4, ":",min:2, " S "...
  139.     end if
  140.  
  141.     if lon > 0.0 then
  142.         deg := floor( lon )
  143.         min := floor( 60 * ( lon - deg ) )
  144.         put deg:4, ":",min:2, " E"
  145.     else
  146.         deg := floor( -lon )
  147.         min := floor( 60 * ( -lon - deg ) )
  148.         put deg:4, ":",min:2, " W"
  149.     end if
  150.     
  151. end procedure
  152.  
  153. %
  154. % 3-point interpolation
  155. %
  156. function interpolate( f0, f1, f2, p : real ) : real
  157.  
  158.     var a, b, f : real
  159.     
  160.     a := f1 - f0
  161.     b := f2 - f1 - a
  162.     f := f0 + p * ( 2 * a + b * ( 2 * p - 1 ) )
  163.  
  164.     return f
  165.  
  166. end function
  167.  
  168. %
  169. % LST at 0h zone time
  170. %
  171. function lst( lon, jd, z : real ) : real
  172.  
  173.     var s, t0 : real
  174.  
  175.     s := 24110.5 + 8640184.812999999 * jd / 36525
  176.     s := s + 86636.6 * z + 86400 * lon
  177.     
  178.     s := s / 86400
  179.     s := s - floor( s )
  180.  
  181.     t0 := s * 360 * DR
  182.  
  183.     return t0
  184.     
  185. end function
  186.  
  187. %
  188. % test an hour for an event
  189. %
  190. procedure test( k, zone : int, 
  191.                 t0, lat, plx : real,
  192.                 var ra, dec, v : rvector )
  193.  
  194.     var ha : rvector
  195.     var a, b, c, d, e, s, z : real
  196.     var hr, min, time : real
  197.     var az, hz, nz, dz : real
  198.     var zchar : string
  199.     var zlet : string := "MLKIHGFEDCBAZNOPQRSTUVWXY"
  200.     label test_exit :
  201.  
  202.     if ra[2] < ra[0] then
  203.         ra[2] := ra[2] + 2 * PI
  204.     end if
  205.     
  206.     ha[0] := t0 - ra[0] + k * K1
  207.     ha[2] := t0 - ra[2] + k * K1 + K1
  208.     
  209.     ha[1]  := ( ha[2] + ha[0] ) / 2         % hour angle at half hour
  210.     dec[1] := ( dec[2] + dec[0] ) / 2       % declination at half hour
  211.  
  212.     s := sin( DR * lat )
  213.     c := cos( DR * lat )
  214.     z := cos( DR * ( 90.567 - 41.685 / plx ) )
  215.  
  216.     if k <= 0 then
  217.         v[0] := s * sin( dec[0] ) + c * cos( dec[0] ) * cos( ha[0] ) - z
  218.     end if
  219.  
  220.     v[2] := s * sin( dec[2] ) + c * cos( dec[2] ) * cos( ha[2] ) - z
  221.     
  222.     if sgn( v[0] ) = sgn( v[2] ) then
  223.         goto test_exit
  224.     end if
  225.     
  226.     v[1] := s * sin( dec[1] ) + c * cos( dec[1] ) * cos( ha[1] ) - z
  227.  
  228.     a := 2 * v[2] - 4 * v[1] + 2 * v[0]
  229.     b := 4 * v[1] - 3 * v[0] - v[2]
  230.     d := b * b - 4 * a * v[0]
  231.  
  232.     if d < 0 then
  233.         goto test_exit
  234.     end if
  235.     
  236.     d := sqrt( d )
  237.     
  238.     if v[0] < 0 and v[2] > 0 then
  239.         put "Moonrise at "...
  240.         Moonrise := true
  241.     end if
  242.  
  243.     if v[0] > 0 and v[2] < 0 then
  244.         put "Moonset at  "...
  245.         Moonset := true
  246.     end if
  247.  
  248.     e := ( -b + d ) / ( 2 * a )
  249.  
  250.     if e > 1 or e < 0 then
  251.         e := ( -b - d ) / ( 2 * a )
  252.     end if
  253.  
  254.     time := k + e + 1 / 120       % round off
  255.  
  256.     hr := floor( time )
  257.     min := floor( ( time - hr ) * 60 )
  258.  
  259.     zchar := "    "
  260.     if abs( zone ) <= 12 then
  261.         zchar[1] := zlet[zone+12]
  262.     end if
  263.  
  264.     put hr:2, ":", min:2, zchar...
  265.  
  266.     % azimuth of the moon at the event
  267.     
  268.     hz := ha[0] + e * ( ha[2] - ha[0] )
  269.     nz := -cos( dec[1] ) * sin( hz )
  270.     dz := c * sin( dec[1] ) - s * cos( dec[1] ) * cos( hz )
  271.     az := arctan( nz / dz ) / DR
  272.  
  273.     if dz < 0 then
  274.         az := az + 180
  275.     end if
  276.     
  277.     if az < 0 then
  278.         az := az + 360
  279.     end if
  280.  
  281.     if az > 360 then
  282.         az := az - 360
  283.     end if
  284.  
  285.     put "azimuth: ", az:5:1
  286.  
  287. test_exit :
  288.  
  289.     v[0] := v[2]
  290.  
  291. end procedure
  292.  
  293. % 
  294. % moon's position using fundamental arguments
  295. %
  296. procedure moon( var ra, dec, plx : real, jd : real )
  297.  
  298.     var d, f, g, h, m, n, s, u, v, w : real 
  299.  
  300.     h := 0.606434 + 0.03660110129 * jd
  301.     m := 0.374897 + 0.03629164709 * jd
  302.     f := 0.259091 + 0.0367481952  * jd
  303.     d := 0.827362 + 0.03386319198 * jd
  304.     n := 0.347343 - 0.00014709391 * jd
  305.     g := 0.993126 + 0.0027377785  * jd
  306.  
  307.     h := h - floor( h )
  308.     m := m - floor( m )
  309.     f := f - floor( f )
  310.     d := d - floor( d )
  311.     n := n - floor( n )
  312.     g := g - floor( g )
  313.  
  314.     h := h * 2 * PI
  315.     m := m * 2 * PI
  316.     f := f * 2 * PI
  317.     d := d * 2 * PI
  318.     n := n * 2 * PI
  319.     g := g * 2 * PI
  320.  
  321.     v := 0.39558 * sin( f + n )
  322.     v := v + 0.082   * sin( f )
  323.     v := v + 0.03257 * sin( m - f - n )
  324.     v := v + 0.01092 * sin( m + f + n )
  325.     v := v + 0.00666 * sin( m - f )
  326.     v := v - 0.00644 * sin( m + f - 2 * d + n )
  327.     v := v - 0.00331 * sin( f - 2 * d + n )
  328.     v := v - 0.00304 * sin( f - 2 * d )
  329.     v := v - 0.0024  * sin( m - f - 2 * d - n )
  330.     v := v + 0.00226 * sin( m + f )
  331.     v := v - 0.00108 * sin( m + f - 2 * d )
  332.     v := v - 0.00079 * sin( f - n )
  333.     v := v + 0.00078 * sin( f + 2 * d + n )
  334.     
  335.     u := 1 - 0.10828 * cos( m )
  336.     u := u - 0.0188  * cos( m - 2 * d )
  337.     u := u - 0.01479 * cos( 2 * d )
  338.     u := u + 0.00181 * cos( 2 * m - 2 * d )
  339.     u := u - 0.00147 * cos( 2 * m )
  340.     u := u - 0.00105 * cos( 2 * d - g )
  341.     u := u - 0.00075 * cos( m - 2 * d + g )
  342.     
  343.     w := 0.10478 * sin( m )
  344.     w := w - 0.04105 * sin( 2 * f + 2 * n )
  345.     w := w - 0.0213  * sin( m - 2 * d )
  346.     w := w - 0.01779 * sin( 2 * f + n )
  347.     w := w + 0.01774 * sin( n )
  348.     w := w + 0.00987 * sin( 2 * d )
  349.     w := w - 0.00338 * sin( m - 2 * f - 2 * n )
  350.     w := w - 0.00309 * sin( g )
  351.     w := w - 0.0019  * sin( 2 * f )
  352.     w := w - 0.00144 * sin( m + n )
  353.     w := w - 0.00144 * sin( m - 2 * f - n )
  354.     w := w - 0.00113 * sin( m + 2 * f + 2 * n )
  355.     w := w - 0.00094 * sin( m - 2 * d + g )
  356.     w := w - 0.00092 * sin( 2 * m - 2 * d )
  357.  
  358.     %    compute right ascension, declination, and parallax
  359.     s  := w / sqrt( u - v * v )
  360.     ra := h + arctan( s / sqrt( 1 - s * s ) )
  361.  
  362.     s   := v / sqrt( u )
  363.     dec := arctan( s / sqrt( 1 - s * s ) )
  364.  
  365.     plx := 60.40974 * sqrt( u )
  366.  
  367. end procedure
  368.  
  369. %
  370. % determine Julian day from calendar date
  371. % ref: Jean Meeus, "Astronomical Algorithms", Willmann-Bell, 1991
  372. %
  373. function julian_day : real
  374.  
  375.     var a, b, day, jd : real
  376.     var month, year : int
  377.     var gregorian : boolean
  378.  
  379.     prompt "Year:"
  380.     get year
  381.     prompt "Month:"
  382.     get month
  383.     prompt "Day:"
  384.     get day
  385.  
  386.     date( floor( day ), month, year )
  387.     
  388.     if year < 1583 then
  389.         gregorian := false
  390.     else
  391.         gregorian := true
  392.     end if
  393.     
  394.     if month = 1 or month = 2 then
  395.         year := year - 1
  396.         month := month + 12
  397.     end if
  398.  
  399.     a := floor( year / 100 )
  400.     if gregorian then
  401.         b := 2 - a + floor( a / 4 )
  402.     else
  403.         b := 0.0
  404.     end if
  405.  
  406.     jd := floor( 365.25 * ( year + 4716 ) ) + 
  407.           floor( 30.6001 * ( month + 1 ) ) + 
  408.           day + b - 1524.5
  409.     
  410.     return jd
  411.     
  412. end function
  413.  
  414. %
  415. % display the date
  416. %
  417. procedure date( d, m, y : int )
  418.  
  419.     var day, month : array[12] of string
  420.     var mi, yi, i : int
  421.  
  422.     month[0]  := "January"
  423.     month[1]  := "February"
  424.     month[2]  := "March"
  425.     month[3]  := "April"
  426.     month[4]  := "May"
  427.     month[5]  := "June"
  428.     month[6]  := "July"
  429.     month[7]  := "August"
  430.     month[8]  := "September"
  431.     month[9]  := "October"
  432.     month[10] := "November"
  433.     month[11] := "December"
  434.  
  435.     day[0] := "Monday"
  436.     day[1] := "Tuesday"
  437.     day[2] := "Wednesday"
  438.     day[3] := "Thursday"
  439.     day[4] := "Friday"
  440.     day[5] := "Saturday"
  441.     day[6] := "Sunday"
  442.  
  443.     mi := m
  444.     yi := y
  445.  
  446.     if y < 1752 then
  447.  
  448.         put month[mi-1], " ", d, ", ", yi
  449.  
  450.     else
  451.  
  452.         if m = 1 or m = 2 then
  453.  
  454.             m := m + 12
  455.             y := y - 1
  456.  
  457.         end if
  458.  
  459.         i := d + 2*m + 3*(m+1) div 5 + y + y div 4 - y div 100 + y div 400
  460.         i := i mod 7
  461.  
  462.         put day[i], ", ", month[mi-1], " ", d, ", ", yi
  463.  
  464.     end if
  465.  
  466. end procedure
  467.  
  468. %
  469. % returns value for sign of argument
  470. %
  471. function sgn( x : real ) : int
  472.  
  473.     var rv : int
  474.     
  475.     if x > 0.0 then
  476.         rv :=  1
  477.     elsif x < 0.0 then
  478.         rv := -1
  479.     else
  480.         rv :=  0
  481.     end if
  482.  
  483.     return rv
  484.  
  485. end function
  486.  
  487. %
  488. % absolute value of argument
  489. %
  490. function abs( x : real ) : real
  491.  
  492.     if x < 0.0 then
  493.         x := -x
  494.     end if
  495.  
  496.     return x
  497.  
  498. end function